Nutrient and chlorophyll analysis
Rosemary Hartman
2022-10-07
Nutrients and chlorophyll
This is the nutrient and chlorophyll analysis for the Drought paper. It pulls the data that Dave organized in the Drought Data package.
First I replaced values below the reporting limit with draws from a uniform distribution, then I averaged by season, region, and year. Then I log-transformed. Sam said he thought log-transforming the mean was better than taking the mean of log-transformed data, and I trust Sam more than me.
Amonia
I’m going to start with Ammonia. My initial look at Ammonia made it seem pretty boring. I’d lean towards leaving it out of the paper, but I”ll see if including some sort of dummy variable for WWTP upgrades helps at all.
Cloern 2019 discussed a lot of Ammonia patterns: Ammonium While the silicate source to San Francisco Bay is river inflow, the ammonium source is primarily municipal wastewater (Jassby 2008). Comparison of silicate and ammonium variability provides an opportunity to explore the differing patterns of nutrients delivered by runoff and point sources. The trends of increasing ammonium were significant for most months (Fig. 6B), so I explored models to explain variability of annual mean ammonium concentration. The best-fitting simple model explained 59% of ammonium variability with two processes— Outflow and wastewater inputs of TKN (Fig. 8D). The size effects of Outflow (−1.0) and TKN loading (+0.8) were comparable (Table 4). The relationship with Outflow is nonlinear and negative, with highest ammonium concentrations during dry years and lowest during wet years (Fig. 8D). This is opposite the Outflow-silicate relationship because river inflow is a source of silicate but it dilutes wastewater-derived ammonium. The shape of nutrient-Outflow relationships can therefore give information about the relative importance of riverine and point-source inputs. Ammonium concentration in the estuary was significantly and linearly related to TKN loading (Fig. 8D), so the long-term trend of increasing ammonium concentration (Fig. 4; Table 2) can be attributed to increasing wastewater loading over time (Fig. 7E), but not to changes in Outflow and dilution.
Ammonium at the monthly scale Analyses of water-quality time series data can yield different insights into processes depending on the time scale considered. As an example, I explored regression models to identify the key drivers of ammonium variability at the monthly scale to compare against analyses at the annual scale. This was motivated by the strong monthly pattern of ammonium trends—largest increases in winter and smallest in summer (Fig. 6B), suggesting a possible influence of seasonal temperature variability. The data series included 346 monthly measurements beginning January 1987. The best fitting model described monthly ammonium variability as a function of three variables—Outflow and TKN loading (as above), and temperature. Monthly ammonium was a nonlinear and negative function of Outflow and a linear positive function of TKN loading (Fig. 9)—the same functional responses observed at the annual scale (Fig. 8D). The Outflow effect was −2.1 and the TKN-loading effect was +1.0. However, at the monthly time scale, temperature had the largest effect (−5.8) through a negative linear relationship. Thus, both annual and monthly variability of ammonium, including trends of increase over time, are associated with variability of TKN loading as a source and Outflow through its dilution effect. However, a stronger temperature effect was revealed by analyzing monthly data. This effect is not evident in analyses of the annual data because annual temperature fluctuations are small relative to the high-amplitude monthly variability (Fig. 4A). The seasonal temperature effect on ammonium concentration in the estuary is presumably a response to the near-linear relationship between nitrification rate and temperature (e.g., Sudarno et al. 2011). Highest nitrification rates during the warm season would explain the negative relationship between ammonium and temperature (Fig. 9A), the seasonal pattern of lowest ammonium concentrations during summer (Fig. 4B), and the seasonal pattern of ammonium increase over time—largest in winter and near-zero in summer (Fig. 6B). Thus, the long-term trends of ammonium increase in the estuary can be attributed to increased wastewater inputs. But the expression of increased loading is not evident in summer (Fig. 6B), implying that wastewater ammonium is nitrified in the time of transit from the wastewater source to the sampling site in the estuary during the warm, low-flow season. Direct measurements of nitrate production and travel time in the lower Sacramento River support this conclusion (Kraus et al. 2017).
#replace values below the reporting limit
Ammonia <- raw_nutr_1975_2021 %>%
# Remove NA values in DissAmmonia
tidyr::drop_na(DissAmmonia) %>%
# Replace <RL values with random value
drt_replace_rl(DissAmmonia, DissAmmonia_Sign) %>%
# Calculate seasonal-regional averages
drt_avg_data(DissAmmonia, avg_type = "both", month_na = "relaxed") %>%
# Add year assignments
drt_add_yr_assign() %>%
mutate(YearType = factor(YearType, levels = c("Critical", "Dry", "Below Normal",
"Above Normal", "Wet")),
LogAm = log(DissAmmonia))
## `geom_smooth()` using formula 'y ~ x'
Now we’ll do the statistics. I’m going to take out 2021, because of the Sac WWTP upgrade. Stockton WWTP was upgraded in 2006
Ammonia = filter(Ammonia, YearAdj < 2021) %>%
mutate(StocktonUpgrade = case_when(YearAdj < 2006 ~ "Before",
YearAdj >= 2006 ~ "After"))
Am1 = lm(LogAm ~ Drought*Region+ Drought*Season + Region*Season, data = Ammonia)
summary(Am1)##
## Call:
## lm(formula = LogAm ~ Drought * Region + Drought * Season + Region *
## Season, data = Ammonia)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.91255 -0.24269 0.03077 0.27067 1.47669
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.56021 0.08046 -31.818 < 2e-16 ***
## DroughtN -0.01760 0.10151 -0.173 0.862372
## DroughtW -0.11432 0.10612 -1.077 0.281741
## RegionNorth 1.26973 0.10679 11.890 < 2e-16 ***
## RegionSouthCentral -0.25081 0.10679 -2.349 0.019137 *
## RegionSuisun Bay -0.08612 0.10454 -0.824 0.410343
## SeasonSpring 0.27482 0.10486 2.621 0.008972 **
## SeasonSummer -0.35780 0.10514 -3.403 0.000706 ***
## SeasonWinter 0.68317 0.10636 6.423 2.54e-10 ***
## DroughtN:RegionNorth -0.14117 0.11003 -1.283 0.199929
## DroughtW:RegionNorth -0.40201 0.11524 -3.488 0.000518 ***
## DroughtN:RegionSouthCentral 0.24124 0.11003 2.193 0.028689 *
## DroughtW:RegionSouthCentral 0.67790 0.11524 5.883 6.40e-09 ***
## DroughtN:RegionSuisun Bay -0.03632 0.10823 -0.336 0.737274
## DroughtW:RegionSuisun Bay 0.04793 0.11318 0.423 0.672086
## DroughtN:SeasonSpring -0.37314 0.10968 -3.402 0.000709 ***
## DroughtW:SeasonSpring -0.51036 0.11446 -4.459 9.67e-06 ***
## DroughtN:SeasonSummer -0.07444 0.10968 -0.679 0.497565
## DroughtW:SeasonSummer 0.06216 0.11529 0.539 0.589968
## DroughtN:SeasonWinter -0.11333 0.11000 -1.030 0.303226
## DroughtW:SeasonWinter -0.35483 0.11507 -3.084 0.002130 **
## RegionNorth:SeasonSpring -0.63148 0.13113 -4.816 1.82e-06 ***
## RegionSouthCentral:SeasonSpring 0.04935 0.13113 0.376 0.706784
## RegionSuisun Bay:SeasonSpring -0.02140 0.12848 -0.167 0.867785
## RegionNorth:SeasonSummer -0.07840 0.13149 -0.596 0.551218
## RegionSouthCentral:SeasonSummer -0.24651 0.13149 -1.875 0.061269 .
## RegionSuisun Bay:SeasonSummer -0.02062 0.12921 -0.160 0.873280
## RegionNorth:SeasonWinter -1.20297 0.13149 -9.149 < 2e-16 ***
## RegionSouthCentral:SeasonWinter 0.40063 0.13149 3.047 0.002404 **
## RegionSuisun Bay:SeasonWinter 0.09303 0.12959 0.718 0.473107
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4333 on 665 degrees of freedom
## (41 observations deleted due to missingness)
## Multiple R-squared: 0.6067, Adjusted R-squared: 0.5895
## F-statistic: 35.37 on 29 and 665 DF, p-value: < 2.2e-16plot(Am1)
visreg(Am1, xvar = "Drought", by = "Region")
visreg(Am1, xvar = "Drought", by = "Season")
visreg(Am1, xvar = "Region", by = "Season")
amp = emmeans(Am1, pairwise ~ Drought*Region)
plot(amp, comparisons = T)## I bet you wanted to call this with just object[[1]] - use '[[]]' or which' if I'm wrong.
## See '? emm_list' for more information
ampx = emmeans(Am1, pairwise ~ Drought*Season)
plot(ampx, comparisons = T)## I bet you wanted to call this with just object[[1]] - use '[[]]' or which' if I'm wrong.
## See '? emm_list' for more information
Now add the treatment plant term
Am2 = lm(LogAm ~ Drought*Region+ Drought*Season + Region*Season + Region*StocktonUpgrade, data = Ammonia)
summary(Am2)##
## Call:
## lm(formula = LogAm ~ Drought * Region + Drought * Season + Region *
## Season + Region * StocktonUpgrade, data = Ammonia)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.77961 -0.22370 0.02605 0.24962 1.04626
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.427775 0.078997 -30.733 < 2e-16
## DroughtN 0.038573 0.091773 0.420 0.674400
## DroughtW -0.001974 0.098839 -0.020 0.984073
## RegionNorth 1.006119 0.105254 9.559 < 2e-16
## RegionSouthCentral -0.728322 0.105254 -6.920 1.07e-11
## RegionSuisun Bay -0.138841 0.104039 -1.335 0.182495
## SeasonSpring 0.273287 0.093650 2.918 0.003641
## SeasonSummer -0.359631 0.093898 -3.830 0.000140
## SeasonWinter 0.673268 0.095010 7.086 3.54e-12
## StocktonUpgradeBefore -0.265380 0.065987 -4.022 6.44e-05
## DroughtN:RegionNorth -0.272378 0.100718 -2.704 0.007020
## DroughtW:RegionNorth -0.653626 0.111290 -5.873 6.77e-09
## DroughtN:RegionSouthCentral -0.009264 0.100718 -0.092 0.926745
## DroughtW:RegionSouthCentral 0.202652 0.111290 1.821 0.069068
## DroughtN:RegionSuisun Bay -0.057905 0.098805 -0.586 0.558041
## DroughtW:RegionSuisun Bay 0.003604 0.108681 0.033 0.973553
## DroughtN:SeasonSpring -0.373724 0.097952 -3.815 0.000149
## DroughtW:SeasonSpring -0.510362 0.102223 -4.993 7.63e-07
## DroughtN:SeasonSummer -0.075031 0.097952 -0.766 0.443953
## DroughtW:SeasonSummer 0.061349 0.102962 0.596 0.551480
## DroughtN:SeasonWinter -0.106400 0.098250 -1.083 0.279224
## DroughtW:SeasonWinter -0.348761 0.102777 -3.393 0.000732
## RegionNorth:SeasonSpring -0.629764 0.117107 -5.378 1.05e-07
## RegionSouthCentral:SeasonSpring 0.051066 0.117107 0.436 0.662936
## RegionSuisun Bay:SeasonSpring -0.018650 0.114747 -0.163 0.870936
## RegionNorth:SeasonSummer -0.076169 0.117434 -0.649 0.516813
## RegionSouthCentral:SeasonSummer -0.244278 0.117434 -2.080 0.037898
## RegionSuisun Bay:SeasonSummer -0.017656 0.115400 -0.153 0.878445
## RegionNorth:SeasonWinter -1.196797 0.117441 -10.191 < 2e-16
## RegionSouthCentral:SeasonWinter 0.406797 0.117441 3.464 0.000567
## RegionSuisun Bay:SeasonWinter 0.096418 0.115752 0.833 0.405161
## RegionNorth:StocktonUpgradeBefore 0.562288 0.094669 5.939 4.62e-09
## RegionSouthCentral:StocktonUpgradeBefore 1.043575 0.094669 11.023 < 2e-16
## RegionSuisun Bay:StocktonUpgradeBefore 0.103600 0.093252 1.111 0.266988
##
## (Intercept) ***
## DroughtN
## DroughtW
## RegionNorth ***
## RegionSouthCentral ***
## RegionSuisun Bay
## SeasonSpring **
## SeasonSummer ***
## SeasonWinter ***
## StocktonUpgradeBefore ***
## DroughtN:RegionNorth **
## DroughtW:RegionNorth ***
## DroughtN:RegionSouthCentral
## DroughtW:RegionSouthCentral .
## DroughtN:RegionSuisun Bay
## DroughtW:RegionSuisun Bay
## DroughtN:SeasonSpring ***
## DroughtW:SeasonSpring ***
## DroughtN:SeasonSummer
## DroughtW:SeasonSummer
## DroughtN:SeasonWinter
## DroughtW:SeasonWinter ***
## RegionNorth:SeasonSpring ***
## RegionSouthCentral:SeasonSpring
## RegionSuisun Bay:SeasonSpring
## RegionNorth:SeasonSummer
## RegionSouthCentral:SeasonSummer *
## RegionSuisun Bay:SeasonSummer
## RegionNorth:SeasonWinter ***
## RegionSouthCentral:SeasonWinter ***
## RegionSuisun Bay:SeasonWinter
## RegionNorth:StocktonUpgradeBefore ***
## RegionSouthCentral:StocktonUpgradeBefore ***
## RegionSuisun Bay:StocktonUpgradeBefore
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3869 on 661 degrees of freedom
## (41 observations deleted due to missingness)
## Multiple R-squared: 0.6882, Adjusted R-squared: 0.6726
## F-statistic: 44.2 on 33 and 661 DF, p-value: < 2.2e-16plot(Am2)
visreg(Am2, xvar = "Drought", by = "Region")
visreg(Am2, xvar = "Drought", by = "Season")
visreg(Am2, xvar = "Region", by = "Season")
visreg(Am2, xvar = "Region", by = "StocktonUpgrade")
amp2 = emmeans(Am2, pairwise ~ Drought*Region)
plot(amp2, comparisons = T)## I bet you wanted to call this with just object[[1]] - use '[[]]' or which' if I'm wrong.
## See '? emm_list' for more information
amp2x = emmeans(Am2, pairwise ~ Drought*Season)
plot(amp2x, comparisons = T)## I bet you wanted to call this with just object[[1]] - use '[[]]' or which' if I'm wrong.
## See '? emm_list' for more information
AIC(Am2)## [1] 687.6603AIC(Am1)## [1] 841.0231Well, the treatment plant term improves the model a bit, but it doesn’t actually change hte drought story. However includeing the region*season interaction and taking out 2021 does make this a little more interesting than before. Droughts in the north are definitely higher in Ammonia.
Now do nitrate
#replace values below the reporting limit
Nitrate <- raw_nutr_1975_2021 %>%
# Remove NA values in DissNitrate
tidyr::drop_na(DissNitrateNitrite) %>%
# Replace <RL values with random value
drt_replace_rl(DissNitrateNitrite, DissNitrateNitrite_Sign) %>%
# Calculate seasonal-regional averages
drt_avg_data(DissNitrateNitrite, avg_type = "both", month_na = "relaxed") %>%
# Add year assignments
drt_add_yr_assign() %>%
mutate(YearType = factor(YearType, levels = c("Critical", "Dry", "Below Normal",
"Above Normal", "Wet")),
LogNat = log(DissNitrateNitrite))Exploritory plots
## `geom_smooth()` using formula 'y ~ x'
OK, now for the statistics. The treatment plant upgrades didn’t have as big an impact on nitrate, so I’ll leave it out for now, but I can put it back in if folks think it’s important.
Nat2 = lm(LogNat ~Drought*Region + Drought*Season + Region*Season, data = Nitrate)
plot(Nat2)
visreg(Nat2, xvar = "Drought", by = "Region")
visreg(Nat2, xvar = "Drought", by = "Season")
visreg(Nat2, xvar = "Region", by = "Season")
summary(Nat2)##
## Call:
## lm(formula = LogNat ~ Drought * Region + Drought * Season + Region *
## Season, data = Nitrate)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.35923 -0.18500 0.00404 0.18472 0.79193
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.017778 0.052548 -19.369 < 2e-16 ***
## DroughtN -0.320539 0.066241 -4.839 1.60e-06 ***
## DroughtW -0.377771 0.069503 -5.435 7.50e-08 ***
## RegionNorth -0.970893 0.068503 -14.173 < 2e-16 ***
## RegionSouthCentral 0.517744 0.068503 7.558 1.25e-13 ***
## RegionSuisun Bay 0.091188 0.068512 1.331 0.183618
## SeasonSpring 0.080514 0.068457 1.176 0.239937
## SeasonSummer -0.167968 0.068457 -2.454 0.014379 *
## SeasonWinter 0.233185 0.069016 3.379 0.000768 ***
## DroughtN:RegionNorth 0.230551 0.070840 3.255 0.001189 **
## DroughtW:RegionNorth 0.303034 0.074546 4.065 5.33e-05 ***
## DroughtN:RegionSouthCentral 0.066926 0.070840 0.945 0.345101
## DroughtW:RegionSouthCentral 0.233552 0.074546 3.133 0.001801 **
## DroughtN:RegionSuisun Bay -0.071598 0.070926 -1.009 0.313088
## DroughtW:RegionSuisun Bay -0.009072 0.074371 -0.122 0.902950
## DroughtN:SeasonSpring -0.021817 0.070753 -0.308 0.757904
## DroughtW:SeasonSpring -0.149871 0.074206 -2.020 0.043790 *
## DroughtN:SeasonSummer -0.023023 0.070753 -0.325 0.744976
## DroughtW:SeasonSummer 0.003852 0.074206 0.052 0.958618
## DroughtN:SeasonWinter 0.252282 0.070928 3.557 0.000400 ***
## DroughtW:SeasonWinter 0.146839 0.074898 1.961 0.050321 .
## RegionNorth:SeasonSpring 0.032742 0.084602 0.387 0.698858
## RegionSouthCentral:SeasonSpring 0.104571 0.084602 1.236 0.216847
## RegionSuisun Bay:SeasonSpring 0.048311 0.084602 0.571 0.568151
## RegionNorth:SeasonSummer -0.281512 0.084602 -3.328 0.000921 ***
## RegionSouthCentral:SeasonSummer -0.176613 0.084602 -2.088 0.037188 *
## RegionSuisun Bay:SeasonSummer 0.105284 0.084602 1.244 0.213734
## RegionNorth:SeasonWinter 0.165879 0.085071 1.950 0.051580 .
## RegionSouthCentral:SeasonWinter 0.219130 0.085071 2.576 0.010199 *
## RegionSuisun Bay:SeasonWinter -0.066024 0.085063 -0.776 0.437902
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.29 on 718 degrees of freedom
## (4 observations deleted due to missingness)
## Multiple R-squared: 0.8208, Adjusted R-squared: 0.8135
## F-statistic: 113.4 on 29 and 718 DF, p-value: < 2.2e-16pairsNN = emmeans(Nat2, pairwise ~ Drought*Region)
plot(pairsNN, comparisons = T)## I bet you wanted to call this with just object[[1]] - use '[[]]' or which' if I'm wrong.
## See '? emm_list' for more information
pairsNN2 = emmeans(Nat2, pairwise ~ Drought*Season)
plot(pairsNN2, comparisons = T)## I bet you wanted to call this with just object[[1]] - use '[[]]' or which' if I'm wrong.
## See '? emm_list' for more information
OK, there is a big main effect of drought here, with droughts having higher nitrate, especially in the spring, summer and fall (no effect in the winter). If I’m reading this right, droughts mean increased nitrate in Suisun Bay and the Confluence, but not the Norht, whereas Ammonia increased in the North. Interesting.
Ortho-phosphate
According to Cloern, total photphorus loading from the WWTP decreased a lot in the 80s, but there was no real phosphorus trend over time.
Phosphate Although phosphate concentration did not change over the long term (Table 2), there was substantial variability of phosphate concentration as an early period of increase during the 1980s, peak concentrations in the early 1990s (Fig. 3) followed by a period of significant decrease (Fig. 5B). Records of phosphorus loading from the Sacramento Regional Wastewater Treatment Facility are not available during the early era of phosphate increase. However, measurements since 1987 showed a steep decrease in wastewater TP loading in the early 1990s that paralleled the downstream phosphate decrease at Sta. D8. This suggests that wastewater effluent is a significant source of phosphorus, as well as ammonium nitrogen, to the upper estuary. The best-fitting regression model (Table 4) explained 73% of annual-mean phosphate concentration as a negative function of Outflow (effect size = −0.74) and a positive function of wastewater TP loading (effect size = +0.25), similar to model results for ammonium. Relative to theOutflow effect, the effect of wastewater TP loading was smaller than the effect of ammonium loading (Table 4). This suggests that the wastewater enrichment of riverine phosphate is smaller than its enrichment of riverine ammonium. The parallel steep declines of phosphorus loading and concentrations in the estuary during the 1990s (Fig. 3) occurred after it became evident that phosphorus enrichment was degrading water quality of U.S. lakes and rivers. State bans on phosphorus-containing detergents grew in the 1970s and 1980s, and by 1994 the manufacture of phosphorus-based detergents ended (Litke 1999). Thus, the oscillating patterns of phosphate concentrations in San Francisco Bay (and other estuaries and lakes) include eras of increase associated with population growth (Jassby 2008) and decrease after policies were implemented to remove phosphorus from household detergents.
Phos <- raw_nutr_1975_2021 %>%
# Remove NA values in DissNitrate
tidyr::drop_na(DissOrthophos) %>%
# Replace <RL values with random value
drt_replace_rl(DissOrthophos, DissOrthophos_Sign) %>%
# Calculate seasonal-regional averages
drt_avg_data(DissOrthophos, avg_type = "both", month_na = "relaxed") %>%
# Add year assignments
drt_add_yr_assign() %>%
mutate(YearType = factor(YearType, levels = c("Critical", "Dry", "Below Normal",
"Above Normal", "Wet")),
LogPhos = log(DissOrthophos)) plots
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
Huh. Cloern didn’t see any overal trends in Phosphorus, but I hink
that’s a pattern. Anyway.
Log-transforming phosphorus didn’t look quite as nice as nitrogen and ammonium. Square-root transformation looked better, but the residuals acutlaly looked pretty good for both, so I just stuck with the log-transformation to keep it consistent.
Phos2 = lm(LogPhos ~Drought*Region + Drought*Season + Region*Season, data = Phos)
plot(Phos2)
visreg(Phos2, xvar = "Drought", by = "Region")
visreg(Phos2, xvar = "Drought", by = "Season")
visreg(Phos2, xvar = "Region", by = "Season")
summary(Phos2)##
## Call:
## lm(formula = LogPhos ~ Drought * Region + Drought * Season +
## Region * Season, data = Phos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.93271 -0.20412 -0.00449 0.20644 0.92730
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.449115 0.055383 -44.221 < 2e-16 ***
## DroughtN -0.178141 0.069819 -2.551 0.01093 *
## DroughtW -0.448946 0.073221 -6.131 1.44e-09 ***
## RegionNorth -0.134761 0.072184 -1.867 0.06232 .
## RegionSouthCentral 0.039391 0.072184 0.546 0.58544
## RegionSuisun Bay 0.179507 0.072213 2.486 0.01315 *
## SeasonSpring -0.230832 0.072155 -3.199 0.00144 **
## SeasonSummer -0.123505 0.072155 -1.712 0.08739 .
## SeasonWinter -0.194480 0.072703 -2.675 0.00764 **
## DroughtN:RegionNorth 0.073187 0.074666 0.980 0.32732
## DroughtW:RegionNorth -0.056897 0.078302 -0.727 0.46768
## DroughtN:RegionSouthCentral 0.094017 0.074666 1.259 0.20838
## DroughtW:RegionSouthCentral 0.200607 0.078302 2.562 0.01061 *
## DroughtN:RegionSuisun Bay -0.021768 0.074757 -0.291 0.77099
## DroughtW:RegionSuisun Bay 0.036316 0.078388 0.463 0.64330
## DroughtN:SeasonSpring -0.100304 0.074575 -1.345 0.17904
## DroughtW:SeasonSpring -0.080454 0.078214 -1.029 0.30400
## DroughtN:SeasonSummer -0.025133 0.074575 -0.337 0.73621
## DroughtW:SeasonSummer 0.003306 0.078214 0.042 0.96630
## DroughtN:SeasonWinter 0.167743 0.074759 2.244 0.02515 *
## DroughtW:SeasonWinter 0.147396 0.078390 1.880 0.06047 .
## RegionNorth:SeasonSpring -0.205459 0.089172 -2.304 0.02150 *
## RegionSouthCentral:SeasonSpring 0.374907 0.089172 4.204 2.95e-05 ***
## RegionSuisun Bay:SeasonSpring -0.059016 0.089172 -0.662 0.50829
## RegionNorth:SeasonSummer -0.188936 0.089172 -2.119 0.03445 *
## RegionSouthCentral:SeasonSummer 0.173622 0.089172 1.947 0.05192 .
## RegionSuisun Bay:SeasonSummer 0.021605 0.089172 0.242 0.80863
## RegionNorth:SeasonWinter -0.162972 0.089416 -1.823 0.06878 .
## RegionSouthCentral:SeasonWinter 0.460494 0.089416 5.150 3.37e-07 ***
## RegionSuisun Bay:SeasonWinter -0.112532 0.089658 -1.255 0.20984
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3057 on 720 degrees of freedom
## (2 observations deleted due to missingness)
## Multiple R-squared: 0.5296, Adjusted R-squared: 0.5107
## F-statistic: 27.96 on 29 and 720 DF, p-value: < 2.2e-16pairsP = emmeans(Phos2, pairwise ~ Drought*Region)
plot(pairsP, comparisons = T)## I bet you wanted to call this with just object[[1]] - use '[[]]' or which' if I'm wrong.
## See '? emm_list' for more information
plot(emmeans(Phos2, pairwise ~ Drought), comparisons = T)## NOTE: Results may be misleading due to involvement in interactions
## I bet you wanted to call this with just object[[1]] - use '[[]]' or which' if I'm wrong.
## See '? emm_list' for more information
For phosphorus, the story is still really the main effect of drought. Lower Phosphorus in wet years.
#now again as a linear model versus sac valley index
Phos2 = lm(LogPhos ~ SVIndex*Region + SVIndex*Season + Region*Season, data = Phos)
plot(Phos2)
visreg(Phos2, xvar = "SVIndex", by = "Region")
visreg(Phos2, xvar = "SVIndex", by = "Season")
visreg(Phos2, xvar = "Region", by = "Season")
summary(Phos2)##
## Call:
## lm(formula = LogPhos ~ SVIndex * Region + SVIndex * Season +
## Region * Season, data = Phos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.04466 -0.18677 -0.00155 0.18192 0.81226
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.065074 0.081180 -25.438 < 2e-16 ***
## SVIndex -0.070965 0.009020 -7.867 1.31e-14 ***
## RegionNorth -0.025468 0.094724 -0.269 0.78811
## RegionSouthCentral -0.124749 0.094724 -1.317 0.18826
## RegionSuisun Bay 0.155678 0.094865 1.641 0.10122
## SeasonSpring -0.165254 0.094581 -1.747 0.08102 .
## SeasonSummer -0.076632 0.094581 -0.810 0.41808
## SeasonWinter -0.315187 0.095399 -3.304 0.00100 **
## SVIndex:RegionNorth -0.013122 0.009650 -1.360 0.17431
## SVIndex:RegionSouthCentral 0.031302 0.009650 3.244 0.00123 **
## SVIndex:RegionSuisun Bay 0.003424 0.009673 0.354 0.72348
## SVIndex:SeasonSpring -0.014920 0.009627 -1.550 0.12161
## SVIndex:SeasonSummer -0.006883 0.009627 -0.715 0.47482
## SVIndex:SeasonWinter 0.027057 0.009673 2.797 0.00529 **
## RegionNorth:SeasonSpring -0.205459 0.081780 -2.512 0.01221 *
## RegionSouthCentral:SeasonSpring 0.374907 0.081780 4.584 5.36e-06 ***
## RegionSuisun Bay:SeasonSpring -0.059016 0.081780 -0.722 0.47074
## RegionNorth:SeasonSummer -0.188936 0.081780 -2.310 0.02115 *
## RegionSouthCentral:SeasonSummer 0.173622 0.081780 2.123 0.03409 *
## RegionSuisun Bay:SeasonSummer 0.021605 0.081780 0.264 0.79171
## RegionNorth:SeasonWinter -0.165025 0.082006 -2.012 0.04455 *
## RegionSouthCentral:SeasonWinter 0.458442 0.082006 5.590 3.21e-08 ***
## RegionSuisun Bay:SeasonWinter -0.112768 0.082227 -1.371 0.17067
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2803 on 727 degrees of freedom
## (2 observations deleted due to missingness)
## Multiple R-squared: 0.6005, Adjusted R-squared: 0.5884
## F-statistic: 49.68 on 22 and 727 DF, p-value: < 2.2e-16Chlorophyll
This is what Cloern has to say:
Chl a Phytoplankton biomass in estuaries is regulated by riverdriven transport processes and the balance between growth and grazing rates. Chl a at Sta. D8 decreased abruptly in 1987 (Table 3), and the largest declines occurred in summer—May through September (Fig. 6C). The best-fitting model of summer Chl a had an interactive effect between Outflow and clam abundance. This interaction was expected because the relationship between Chl a and Outflow changed after the Potamocorbula introduction (Jassby 2008). Most (80%) of that variability was associated with Outflow and clam abundance treated as a categorical variable (Fig. 8B). No other variables, such as temperature, turbidity (Secchi depth), nutrient concentrations, or loadings had comparable effects on summer Chl a. Summer phytoplankton blooms occurred regularly in this region of the estuary during the 1970s and 1980s. These were attributed to the accumulation of phytoplankton cells withinan estuarine turbidity maximum that was positioned in Suisun Bay (Sta. D8) when Outflow was in the range 100–350 m3 s−1 (Cloern et al. 1983). Thus, flow-phytoplankton relationships were established early in the study of San Francisco Bay. Observations made over the subsequent 3+ decades confirm that summer Chl a decreases when Outflow exceeds about200–300 m3 s−1 (Fig. 8B), primarily through its effect on residence time (Jassby 2008). A summer bloom did not develop in 1977, a year of extreme drought when salinity increased high enough that the upper estuary was colonized by marine benthic filter feeders including M. arenaria (Nichols 1985). The 1977 anomaly was the first evidence that bivalve grazing can suppress blooms in this estuary. Summer blooms disappeared completely after P. amurensis became established as a permanent resident in 1987. The regression model (Table 4) shows that the effect of clam presence was large (−13.1) compared to the outflow effect (−2.5). The nearly fourfold loss of phytoplankton biomass (Table 2) and production restructured biological communities and altered nutrient dynamics, establishing what has become an iconic example of abrupt disturbance by the introduction of a non-native species (Cloern and Jassby 2012). Regression models allow us to compare the effects of that human disturbance with the effect of varying freshwater inflow.
#replace values below the reporting limit
ChlaA <- raw_chla_1975_2021 %>%
# Remove NA values in DissAmmonia
tidyr::drop_na(Chlorophyll) %>%
# Replace <RL values with random value
drt_replace_rl(Chlorophyll, Chlorophyll_Sign) %>%
# Calculate seasonal-regional averages
drt_avg_data(Chlorophyll, avg_type = "both", month_na = "relaxed") %>%
# Add year assignments
drt_add_yr_assign() %>%
mutate(YearType = factor(YearType, levels = c("Critical", "Dry", "Below Normal",
"Above Normal", "Wet")),
LogChl = log(Chlorophyll),
Regime = case_when(YearAdj > 1987 ~ "post-clam",
YearAdj <1988 ~ "pre-clam"))Now some plots
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
And now for some statistics
Chl1 = lm(LogChl ~ Drought*Region+ Drought*Season + Region*Season + Region*Regime, data = ChlaA)
summary(Chl1)##
## Call:
## lm(formula = LogChl ~ Drought * Region + Drought * Season + Region *
## Season + Region * Regime, data = ChlaA)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.95368 -0.31348 -0.02377 0.28078 1.93400
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.818998 0.090823 9.018 < 2e-16 ***
## DroughtN -0.008139 0.115973 -0.070 0.944067
## DroughtW 0.012771 0.120054 0.106 0.915310
## RegionNorth -0.730709 0.118558 -6.163 1.19e-09 ***
## RegionSouthCentral 1.203597 0.118558 10.152 < 2e-16 ***
## RegionSuisun Bay -0.310979 0.118613 -2.622 0.008933 **
## SeasonSpring 0.576438 0.117325 4.913 1.11e-06 ***
## SeasonSummer 0.347300 0.117325 2.960 0.003176 **
## SeasonWinter -0.419778 0.117372 -3.576 0.000372 ***
## Regimepre-clam 0.692098 0.084413 8.199 1.12e-15 ***
## DroughtN:RegionNorth 0.147293 0.125964 1.169 0.242663
## DroughtW:RegionNorth 0.224310 0.129194 1.736 0.082953 .
## DroughtN:RegionSouthCentral -0.302435 0.125964 -2.401 0.016606 *
## DroughtW:RegionSouthCentral -0.397518 0.129194 -3.077 0.002171 **
## DroughtN:RegionSuisun Bay 0.300453 0.126081 2.383 0.017431 *
## DroughtW:RegionSuisun Bay 0.257865 0.129316 1.994 0.046523 *
## DroughtN:SeasonSpring -0.173823 0.121259 -1.433 0.152154
## DroughtW:SeasonSpring -0.297770 0.127177 -2.341 0.019485 *
## DroughtN:SeasonSummer 0.071331 0.121259 0.588 0.556544
## DroughtW:SeasonSummer -0.170360 0.127177 -1.340 0.180816
## DroughtN:SeasonWinter -0.398930 0.121408 -3.286 0.001066 **
## DroughtW:SeasonWinter -0.522755 0.127319 -4.106 4.49e-05 ***
## RegionNorth:SeasonSpring 0.321526 0.144993 2.218 0.026900 *
## RegionSouthCentral:SeasonSpring -0.179922 0.144993 -1.241 0.215050
## RegionSuisun Bay:SeasonSpring -0.064102 0.144993 -0.442 0.658550
## RegionNorth:SeasonSummer 0.108072 0.144993 0.745 0.456300
## RegionSouthCentral:SeasonSummer 0.333124 0.144993 2.298 0.021877 *
## RegionSuisun Bay:SeasonSummer 0.057263 0.144993 0.395 0.693010
## RegionNorth:SeasonWinter 1.064912 0.144993 7.345 5.63e-13 ***
## RegionSouthCentral:SeasonWinter -0.045488 0.144993 -0.314 0.753823
## RegionSuisun Bay:SeasonWinter 0.177521 0.145392 1.221 0.222495
## RegionNorth:Regimepre-clam -0.454639 0.119379 -3.808 0.000152 ***
## RegionSouthCentral:Regimepre-clam -0.320410 0.119379 -2.684 0.007443 **
## RegionSuisun Bay:Regimepre-clam 0.245044 0.119397 2.052 0.040498 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.497 on 717 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.6841, Adjusted R-squared: 0.6695
## F-statistic: 47.05 on 33 and 717 DF, p-value: < 2.2e-16plot(Chl1)
library(car)##
## Attaching package: 'car'## The following object is masked from 'package:dplyr':
##
## recode## The following object is masked from 'package:purrr':
##
## someqqPlot(Chl1)
## [1] 238 328hist(residuals(Chl1))
#better
hist(log(ChlaA$Chlorophyll))
visreg(Chl1, xvar = "Drought", by = "Region")
visreg(Chl1, xvar = "Drought", by = "Season")
visreg(Chl1, xvar = "Region", by = "Season")
pairsC = emmeans(Chl1, pairwise ~ Drought*Region)
plot(pairsC, comparisons = T)## I bet you wanted to call this with just object[[1]] - use '[[]]' or which' if I'm wrong.
## See '? emm_list' for more information
pairsC2 = emmeans(Chl1, pairwise ~ Drought*Season)
plot(pairsC2, comparisons = T)## I bet you wanted to call this with just object[[1]] - use '[[]]' or which' if I'm wrong.
## See '? emm_list' for more information
pairsC3 = emmeans(Chl1, pairwise ~Regime*Region)
plot(pairsC3, comparisons = T)## I bet you wanted to call this with just object[[1]] - use '[[]]' or which' if I'm wrong.
## See '? emm_list' for more information
Well, I added the regime shift at the clam invasion and now my diagnostic plots look better, droughts definitely have higher chlorophyll in the south-central region, and definitely have higher winter chlorophyll YAY!!
Try the linear model with Sac valley index
Chl2 = lm(LogChl ~ SVIndex*Region+ SVIndex*Season + Region*Season + Region*Regime, data = ChlaA)
summary(Chl2)##
## Call:
## lm(formula = LogChl ~ SVIndex * Region + SVIndex * Season + Region *
## Season + Region * Regime, data = ChlaA)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.82351 -0.31874 -0.04438 0.28400 1.89314
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.723453 0.145803 4.962 8.71e-07 ***
## SVIndex 0.013673 0.016213 0.843 0.399329
## RegionNorth -0.758480 0.170309 -4.454 9.78e-06 ***
## RegionSouthCentral 1.582441 0.170309 9.292 < 2e-16 ***
## RegionSuisun Bay -0.384305 0.170586 -2.253 0.024567 *
## SeasonSpring 0.783010 0.169551 4.618 4.58e-06 ***
## SeasonSummer 0.398841 0.169551 2.352 0.018923 *
## SeasonWinter -0.217558 0.169809 -1.281 0.200536
## Regimepre-clam 0.656238 0.082534 7.951 7.10e-15 ***
## SVIndex:RegionNorth 0.015256 0.017383 0.878 0.380428
## SVIndex:RegionSouthCentral -0.071978 0.017383 -4.141 3.87e-05 ***
## SVIndex:RegionSuisun Bay 0.027097 0.017421 1.555 0.120283
## SVIndex:SeasonSpring -0.043008 0.017257 -2.492 0.012918 *
## SVIndex:SeasonSummer -0.009491 0.017257 -0.550 0.582505
## SVIndex:SeasonWinter -0.058462 0.017299 -3.379 0.000765 ***
## RegionNorth:SeasonSpring 0.321526 0.146604 2.193 0.028613 *
## RegionSouthCentral:SeasonSpring -0.179922 0.146604 -1.227 0.220122
## RegionSuisun Bay:SeasonSpring -0.064102 0.146604 -0.437 0.662065
## RegionNorth:SeasonSummer 0.108072 0.146604 0.737 0.461258
## RegionSouthCentral:SeasonSummer 0.333124 0.146604 2.272 0.023362 *
## RegionSuisun Bay:SeasonSummer 0.057263 0.146604 0.391 0.696212
## RegionNorth:SeasonWinter 1.064912 0.146604 7.264 9.74e-13 ***
## RegionSouthCentral:SeasonWinter -0.045488 0.146604 -0.310 0.756440
## RegionSuisun Bay:SeasonWinter 0.176784 0.147009 1.203 0.229547
## RegionNorth:Regimepre-clam -0.417387 0.116720 -3.576 0.000372 ***
## RegionSouthCentral:Regimepre-clam -0.359528 0.116720 -3.080 0.002146 **
## RegionSuisun Bay:Regimepre-clam 0.310076 0.116755 2.656 0.008086 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5025 on 724 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.6739, Adjusted R-squared: 0.6622
## F-statistic: 57.54 on 26 and 724 DF, p-value: < 2.2e-16plot(Chl2)
visreg(Chl2, xvar ="SVIndex", by = "Region")
visreg(Chl2, xvar = "SVIndex", by = "Season")
visreg(Chl2, xvar = "Region", by = "Season")
table(raw_chla_1975_2021$Source, raw_chla_1975_2021$YearAdj)
table(raw_nutr_1975_2021$Source, raw_nutr_1975_2021$YearAdj)